
# Copyright (c) 2018-2021, Carnegie Mellon University
# See LICENSE for details


# ************************************************************
#   Hashing for the lifting schemes (uses HashDP)
# ************************************************************
HashLookupWav := HashLookupDefault;
HashIndexWavRec := function( l )
   local index, param, coef;
   index :=13;
   for param in l[1] do
      for coef in param do
         if IsInt(param) then
         index := index * (param+1);
     elif IsRat(param) then
         index := index * Int(coef);
     elif IsFloat(param) then
         index := index * IntFloat(coef);
         elif IsCyc(param) then
            index := index + 1;
         fi;
      od;
   od;
   return index;
end;

HashIndexWav := (el, hsize) -> (HashIndexWavRec(el) mod hsize) + 1;
HashEqualWav := (l1,l2) -> l1=l2;
HashTableWavelets:=HashTable(HashIndexWav, HashEqualWav);

# *************************************************
#   Library of wavelets and lifting schemes
# *************************************************

# Daubechies 9-7

z:=Indeterminate(Cyclotomics);
z.name:="z";
Local.hc := [0.0378287985799179, -0.0238492975158599, -0.110624027489511, 0.377402688109134, 0.852698653219295, 0.377402688109134, -0.110624027489511, -0.0238492975158599, 0.0378287985799179];;

gc := [0.0645389087111023, -0.0406897526165987, -0.418092440725732, 0.788484872206183, -0.418092440725732, -0.0406897526165987,0.0645389087111023];;
gc:=-gc;


Local.h:=Polynomial(Cyclotomics,hc,-4);
g:=Polynomial(Cyclotomics,gc,-2);

h.name:="Daub97lo";
g.name:="Daub97hi";
g.operations := CopyFields(g.operations, rec(Print:=self >> Print(self.name)));
h.operations := CopyFields(h.operations, rec(Print:=self >> Print(self.name)));

Daub97 := [h,g];

#LS:=List([1..4], k->Poly([Random([-100..100]),Random([-100..100])], - (k mod 2)));
#LS:=Concat([0, Poly([Random([-100..100])],0),  Poly([Random([-100..100])],0)],LS);

# Use only one LS (the only linear phase LS)
 LS:=[[ 0, Poly([-0.86986325969200529],0), Poly([1.1496043053581531],0),
  Poly([0.44350482244527795, 0.44350482244527795], -1),
  Poly([0.88293362717898749,0.88293362717898749], 0),
  Poly([-0.052978640032586351,-0.052978640032586351],-1),
  Poly([-1.5861598671726733,-1.5861598671726733],0) ] ];

# or generate all 27 different schemes
#
# M:=Polyphase(h,g);
# LS:=LiftingScheme(M);
 HashAdd(HashTableWavelets, ListPolyMat([h,g]), LS);



# Rational symmetric

hc:=[-1/8, 1/4, 3/4, 1/4, -1/8];
gc:=[1/4, -1/2, 1/4];

 z:=Indeterminate(Rationals);
 z.name:="z";

h:=Polynomial(Rationals,hc,-2);
g:=Polynomial(Rationals,gc,0);
g.operations.Print:=self >> Print(self.name);
h.operations.Print:=self >> Print(self.name);

h.name:="rat53lo";
g.name:="rat53hi";

linwav:=[h,g];
Local.M:=Polyphase(h,g);
LS := LiftingScheme(M);
HashAdd(HashTableWavelets, ListPolyMat([h,g]), LS);


# Haar
 z:=Indeterminate(Cyclotomics);
 z.name:="z";

 h := 1+z;
 g := 1/2 - 1/2*z;
 Haar_n:=[h,g];
M:=Polyphase(h,g);
LS := LiftingScheme(M);
HashAdd(HashTableWavelets, ListPolyMat([h,g]), LS);


 h := 1+z;
 g := 1-z;
 Haar:=[h,g];
M:=Polyphase(h,g);
LS := LiftingScheme(M);
HashAdd(HashTableWavelets, ListPolyMat([h,g]), LS);


# Daubechies D4

 z:=Indeterminate(Cyclotomics);
 z.name:="z";

# h filter
hc:=([(1-Sqrt(3)),(3-Sqrt(3)),(3+Sqrt(3)),(1+Sqrt(3))]/(4*Sqrt(2)));
# g filter
gc:=([(1+Sqrt(3)),-(3+Sqrt(3)),(3-Sqrt(3)),-(1-Sqrt(3))]/(4*Sqrt(2)));

h:=Polynomial(Cyclotomics,hc,0);
g:=Polynomial(Cyclotomics,gc,-2);

Daub4:=[h,g];
M:=Polyphase(h,g);
LS := LiftingScheme(M);
HashAdd(HashTableWavelets, ListPolyMat([h,g]), LS);

# Daubechies D6
#---------------



hc := [
 3.326705529500826159985115891390056300129233992450683597084705e-01,
 8.068915093110925764944936040887134905192973949948236181650920e-01,
 4.598775021184915700951519421476167208081101774314923066433867e-01,
-1.350110200102545886963899066993744805622198452237811919756862e-01,
-8.544127388202666169281916918177331153619763898808662976351748e-02,
 3.522629188570953660274066471551002932775838791743161039893406e-02,
];

gc := Reversed(hc);
for i in [1..Length(gc)] do
  if (i mod 2 = 0) then gc[i]:=-gc[i]; fi;
od;

Daub6 := [hc,gc];

#h:=Polynomial(Floats,hc,-2);
#g:=Polynomial(Floats,gc,-6);



# Daubechies D8
#---------------


hc := [
 2.303778133088965008632911830440708500016152482483092977910968e-01,
 7.148465705529156470899219552739926037076084010993081758450110e-01,
 6.308807679298589078817163383006152202032229226771951174057473e-01,
 -2.798376941685985421141374718007538541198732022449175284003358e-02,
 -1.870348117190930840795706727890814195845441743745800912057770e-01,
 3.084138183556076362721936253495905017031482172003403341821219e-02,
 3.288301166688519973540751354924438866454194113754971259727278e-02,
 -1.059740178506903210488320852402722918109996490637641983484974e-02
];

gc := Reversed(hc);
for i in [1..Length(gc)] do
  if (i mod 2 = 0) then gc[i]:=-gc[i]; fi;
od;

Daub8 := [hc,gc];

#h:=Polynomial(Floats,hc,-2);
#g:=Polynomial(Floats,gc,-6);

# Daubechies D12
#---------------



hc := [
 1.115407433501094636213239172409234390425395919844216759082360e-01,
 4.946238903984530856772041768778555886377863828962743623531834e-01,
 7.511339080210953506789344984397316855802547833382612009730420e-01,
 3.152503517091976290859896548109263966495199235172945244404163e-01,
-2.262646939654398200763145006609034656705401539728969940143487e-01,
-1.297668675672619355622896058765854608452337492235814701599310e-01,
 9.750160558732304910234355253812534233983074749525514279893193e-02,
 2.752286553030572862554083950419321365738758783043454321494202e-02,
-3.158203931748602956507908069984866905747953237314842337511464e-02,
 5.538422011614961392519183980465012206110262773864964295476524e-04,
 4.777257510945510639635975246820707050230501216581434297593254e-03,
-1.077301085308479564852621609587200035235233609334419689818580e-03
];

gc := Reversed(hc);
for i in [1..Length(gc)] do
  if (i mod 2 = 0) then gc[i]:=-gc[i]; fi;
od;

Daub12 := [hc,gc];

#h:=Polynomial(Floats,hc,-6);
#g:=Polynomial(Floats,gc,-6);


# Daubechies D16
#---------------



hc := [
 5.441584224310400995500940520299935503599554294733050397729280e-02,
 3.128715909142999706591623755057177219497319740370229185698712e-01,
 6.756307362972898068078007670471831499869115906336364227766759e-01,
 5.853546836542067127712655200450981944303266678053369055707175e-01,
-1.582910525634930566738054787646630415774471154502826559735335e-02,
-2.840155429615469265162031323741647324684350124871451793599204e-01,
 4.724845739132827703605900098258949861948011288770074644084096e-04,
 1.287474266204784588570292875097083843022601575556488795577000e-01,
-1.736930100180754616961614886809598311413086529488394316977315e-02,
-4.408825393079475150676372323896350189751839190110996472750391e-02,
 1.398102791739828164872293057263345144239559532934347169146368e-02,
 8.746094047405776716382743246475640180402147081140676742686747e-03,
-4.870352993451574310422181557109824016634978512157003764736208e-03,
-3.917403733769470462980803573237762675229350073890493724492694e-04,
 6.754494064505693663695475738792991218489630013558432103617077e-04,
-1.174767841247695337306282316988909444086693950311503927620013e-04
];

gc := Reversed(hc);
for i in [1..Length(gc)] do
  if (i mod 2 = 0) then gc[i]:=-gc[i]; fi;
od;

Daub16 := [hc,gc];
#h:=Polynomial(Floats,hc,-8);
#g:=Polynomial(Floats,gc,-8);




# Daubechies D30
# z:=Indeterminate(Floats);
# z.name:="z";

hc:= [ 6.133360e-08 , -6.316882e-07 , 1.811270e-06 ,
3.362987e-06 , -2.813330e-05 , 2.579270e-05 , 1.558965e-04 ,
-3.595652e-04 , -3.734824e-04 , 1.943324e-03 , -2.417565e-04 ,
-6.487735e-03 , 5.101000e-03 , 1.508392e-02 , -2.081005e-02 ,
-2.576701e-02 , 5.478055e-02 , 3.387714e-02 , -1.111209e-01 ,
-3.966618e-02 , 1.901467e-01 , 6.528295e-02 , -2.888826e-01 ,
-1.932041e-01 , 3.390025e-01 , 6.458131e-01 , 4.926318e-01 ,
2.060239e-01 , 4.674339e-02 , 4.538537e-03  ];

gc:= [ -4.538537e-03 , 4.674339e-02 , -2.060239e-01 ,
4.926318e-01 , -6.458131e-01 , 3.390025e-01 , 1.932041e-01 ,
-2.888826e-01 , -6.528295e-02 , 1.901467e-01 , 3.966618e-02 ,
-1.111209e-01 , -3.387714e-02 , 5.478055e-02 , 2.576701e-02 ,
-2.081005e-02 , -1.508392e-02 , 5.101000e-03 , 6.487735e-03 ,
-2.417565e-04 , -1.943324e-03 , -3.734824e-04 , 3.595652e-04 ,
1.558965e-04 , -2.579270e-05 , -2.813330e-05 , -3.362987e-06 ,
1.811270e-06 , 6.316882e-07 , 6.133360e-08  ];

h:=Polynomial(Cyclotomics,hc,-14);
g:=Polynomial(Cyclotomics,gc,-14);

Daub30 := [h,g];

LS:=List([1..15], k->Poly([Random([-100..100]),Random([-100..100])], - (k mod 2)));
LS:=Concat([0, Poly([Random([-100..100])],0),  Poly([Random([-100..100])],0)],LS);
HashAdd(HashTableWavelets, ListPolyMat([h,g]), [LS]);



#lifting:= [ 0, (0.097094733608324074181)*z^0,
#  ((0.29548526964299320907))*z + ((1.9693135582762033575)),
#  ((-0.19310907414953903949))*z + ((-0.75349025458655549681)),
#  ((6.418850548852764959))*z + ((5.6046209933038362294)),
#  ((2.2768252163984019631))*z + ((2.8590859893363669286)),
#  ((0.42436816698996654429))*z + ((-0.24306346476809917445)),
#  ((32.824621821962480794))*z + ((7.9226055123935470448)),
#  ((0.14922683209603376797))*z + ((-0.056498191876760694985)),
#  ((57.543585676762852188))*z + ((-1.5098853311803024368)),
#  ((-1.8344976647693662652))*z + ((-5.6547891332523283481)),
#  ((-0.052782687371178790836))*z + ((0.17584142116743275985)),
#  ((-3102.2613278850230927))*z + ((1012.1832446864441408)),
#  ((-0.018859084930037840061))*z + ((0.013204658033145515172)),
#  ((-252.58258260810293905))*z + ((-70.090604941421375429)),
#  ((0.2412682849948278585))*z + ((-0.101870749823854162)), ((-4.1447634142093852105))*z^(
#    -1) + ((-1.7500444501621712501))*z^(
#    -2) + ((-0.80693505491745187719))*z^(
#    -3) + ((-0.35055632777493950236))*z^(
#    -4) + ((-0.1442471074569136591))*z^(
#    -5) + ((-0.055639926391646214732))*z^(
#    -6) + ((-0.019909460195237267677))*z^(
#    -7) + ((-0.00651559931905414981))*z^(
#    -8) + ((-0.001913788821165934113))*z^(
#    -9) + ((-0.00049136620953624493384))*z^(
#    -10) + ((-0.00010605600372416396462))*z^(-11),
#  (1.357366318938735191e-06)*z^0,
#  (736720.7997188672889)*z^0 ]
